#############################################################################
##
#W  groupsearch.gi           DifSets Package                     Dylan Peifer
##
##  Functions set up the algorithm and pass work off to the refining and
##  equivalent list functions for each stage.
##

InstallGlobalFunction( ExportMultipliers, function (ordera, orderb, filename)
    local order, idx, G, sizes, size, k, lambda, multipliers, multiplier, tracksOfGroup;

    PrintTo(filename, "order, index, e, lambda, multiplier, tracks\n");

    for order in [ordera..orderb] do
        for idx in [1..NumberSmallGroups(order)] do
            G := SmallGroup(order, idx);;
            if not IsAbelian(G) then
                continue;;
            fi;;
            sizes := PossibleDifferenceSetPairSizes(G);;
            for size in sizes do
                k := size[3];
                lambda := size[4];
                multipliers := Multipliers(G, k, lambda);
                for multiplier in multipliers do
                    tracksOfGroup := TracksOfGroup(G, multiplier);
                    AppendTo(filename,order,", ",idx,", ",k,", ",lambda,", ",multiplier,", ",tracksOfGroup,"\n");
                od;;
            od;;
        od;;
    od;;
end );

#############################################################################
##
#F  RefiningSeries( <G> )
##
InstallGlobalFunction( RefiningSeries, function (G)
    local bestN, bestsize, N;

    # find a nontrivial normal subgroup of smallest possible size
    bestN := G;
    bestsize := Size(G);

    for N in NormalSubgroups(G) do
        if Size(N) < bestsize and Size(N) > 1 then
            bestN := N;
            bestsize := Size(N);
        fi;
    od;

    # return a chief series ending in this minimal normal subgroup
    return ChiefSeriesThrough(G, [bestN]);
end );

#############################################################################
##
#F  PossibleDifferenceSetSizes( <G> )
##
InstallGlobalFunction( PossibleDifferenceSetSizes, function (G)
    local v, ks, IsSquareFree, Test;

    v := Size(G);
    ks := List([2..Int(v/2)]); # ignore size 1 and larger of complements

    # counting test
    ks := Filtered(ks, k->IsInt(k*(k-1)/(v-1)));

    # BRC test
    if IsEvenInt(v) then
        ks := Filtered(ks, k->IsSquareInt(k - k*(k-1)/(v-1)));
    else
        IsSquareFree := x -> Product(Set(FactorsInt(x))) = x;
        Test := function (k)
            local lambda, a, b, d;
            lambda := k*(k-1)/(v-1);
            a := k - lambda;
            b := (-1)^((v-1)/2) * lambda;
            if IsSquareFree(a) and IsSquareFree(b) then
                d := GcdInt(a, b);
                return Legendre(a, AbsInt(b)) = 1 and
                       Legendre(b, AbsInt(a)) = 1 and
                       Legendre(-a*b/d^2, d) = 1;
            else
                return true;
            fi;
            end;

        ks := Filtered(ks, k->Test(k));
    fi;

    return ks;
end );

#############################################################################
##
#F  PossibleDifferenceSetPairSizes( <G> )
##
# FEATURE PossibleDifferenceSetPairSizes(G)
# k1*k2 = λ*(v-1) + e
# 0 < k1 <= Int(v/2)  # only the smaller of each complementary pair of difference set pairs is include.
# 0 < k2 <= Int(v/2)
# 0 <= e <= Min(k1,k2)
# return as [k1, k2, e, lambda]
InstallGlobalFunction( PossibleDifferenceSetPairSizes, function (G)
    local v, k1, k2, e, lambda, min, ks, item, rs;
    ks:=[];
    v:=Size(G);
    for k1 in [1..Int(v)] do
        for k2 in [1..Int(v)] do
            min := k1;  # min = Min(k1,k2);
            if min>k2 then
                min :=k2;
            fi;
            for e in [0..min] do
                # check if lambda is int 
                lambda:= (k1*k2-e)/(v-1);
                if IsInt(lambda) and lambda > 0 and (e-lambda)>1 then
                    # TODO add additional conditions if needed: n = e - lambda > 1
                    # possible k1&k2, save to ks;
                    Add(ks,[k1, k2, e, lambda]);
                fi;
            od;
        od;
    od;

    rs := [];
    # Filter out parameters for k1 and k2 exchanges, such as [ 2, 3, 0, 1 ], [ 3, 2, 0, 1 ]
    while Length(ks) > 0 do
        # remove a size pair and save to rs
        item := Remove(ks);
        Add(rs,item);

        # remove anything equivalent to item from ks
        RemoveSet(ks,[item[2],item[1],item[3],item[4]]);
    od;
    Sort(rs);
    return rs;
end );


#############################################################################
##
#F  PossibleDivisibleDifferenceSetSizes( <G> )
##
# FEATURE PossibleDivisibleDifferenceSetSizes(G)
# l*l = l+lambda1*(v-1)+lambda2*(w*v-v)
# |G| = w*v
# |N| = v
# N = {0, w, 2w, ..., (v-1)w}
# |D| = l
# return as [w, v, l, lambda1, lambda2]
InstallGlobalFunction( PossibleDivisibleDifferenceSetSizes, function (G)
    local n, factors, factorCombs, combs, idx, w, v, l, lambda1, lambda2, ks;

    n:=Size(G);
    factors := Factors(n);
    factorCombs :=[];
    for idx in [1..(Length(factors)-1)] do
        combs := Combinations(factors, idx);
        Apply(combs , x->Product(x));
        Append(factorCombs, combs);
    od;
    factorCombs := Set(factorCombs);  # deduplication

    ks := []; # store possible sizes
    # we iterate the list factorCombs, so that we have w&v
    for w in factorCombs do
        v := n / w;
        for l in [2..Int(n/2)] do  # Only the smaller of each complementary pair of divisible difference sets is included
            for lambda1 in [1..Int(l*(l-1)/(v-1))] do 
                lambda2 := (l*(l-1)-lambda1*(v-1))/(w*v - v);
                if lambda2 > 0 and IsInt(lambda2) then
                    # one possible size
                    Add(ks,[w, v, l, lambda1, lambda2]);
                fi;
            od;
        od;

    od;
    return ks;
end );

#############################################################################
##
#F  PossibleDivisibleDifferenceSetPairSizes( <G> )
##
# FEATURE PossibleDivisibleDifferenceSetPairSizes(G)
# l1*l2 = le+lambda1*(v-1)+lambda2*(w*v-v)
# |G| = w*v
# |N| = v
# N = {0, w, 2w, ..., (v-1)w}
# |D1| = l1, |D2| = l2, |D1∩D2| = le
# return as [w, v, l1, l2, le, lambda1, lambda2]
InstallGlobalFunction( PossibleDivisibleDifferenceSetPairSizes, function (G)
    local n, factors, factorCombs, combs, idx, w, v, l1, l2, min, le, lambda1, lambda2, ks, item, rs;

    n:=Size(G);
    factors := Factors(n);
    factorCombs :=[];
    for idx in [1..(Length(factors)-1)] do
        combs := Combinations(factors, idx);
        Apply(combs , x->Product(x));
        Append(factorCombs, combs);
    od;
    factorCombs := Set(factorCombs);  # deduplication

    ks := []; # store possible sizes
    # we iterate the list factorCombs, so that we have w&v
    for w in factorCombs do
        v := n / w;
        for l1 in [2..Int(n/2)] do  # Only the smaller of each complementary pair of divisible difference sets is included
            for l2 in [2..Int(n/2)] do  # Only the smaller of each complementary pair of divisible difference sets is included
                min := l1;
                if min > l2 then
                    min := l2;
                fi;
                for le in [0..min] do
                    for lambda1 in [1..Int((l1*l2-le)/(v-1))] do 
                        lambda2 := (l1*l2-le-lambda1*(v-1))/(w*v - v);
                        if lambda2 > 0 and IsInt(lambda2) then
                            # one possible size
                            Add(ks,[w, v, l1, l2, le, lambda1, lambda2]);
                        fi;
                    od;
                od;
            od;
        od;

    od;

    # Filter out parameters for k1 and k2 exchanges
    rs := [];
    while Length(ks)>0 do
        item := Remove(ks);
        Add(rs,item);

        RemoveSet(ks,[item[1],item[2],item[4],item[3],item[5],item[6],item[7]]);
    od;
    Sort(rs);
    return rs;
end );

#############################################################################
##
#F  DifferenceSetsOfSizeK( <G>, <k> )
##
InstallGlobalFunction( DifferenceSetsOfSizeK, function (G, k)
    local series, oldN, difsums, i, newN, difsets;

    series := RefiningSeries(G);

    # at start the quotient group G/G is trivial and the only difsum is k
    oldN := G;
    difsums := [[k]];

    # perform search and remove equivalents for each stage
    for i in [2..Length(series)-1] do

        newN := series[i];
        difsums := SomeRefinedDifferenceSums(G, oldN, newN, difsums);
        difsums := EquivalentFreeListOfDifferenceSums(G, newN, difsums);

        oldN := newN; # prepare for next stage

    od;

    # in final stage refine to difsets and remove equivalents
    # difsets := SomeRefinedDifferenceSets(G, oldN, difsums);
    difsets := SomeRefinedDifferenceSets(G, oldN, difsums);
    difsets := SmallestEquivalentFreeListOfDifferenceSets(G, difsets);

    return difsets;
end );


#############################################################################
##
#F  DifferenceSetsOfSizeKOptimized( <G>, <k> )
##
InstallGlobalFunction( DifferenceSetsOfSizeKOptimized, function (G, k)
    local series, oldN, difsums, i, newN, difsets;

    series := RefiningSeries(G);

    # at start the quotient group G/G is trivial and the only difsum is k
    oldN := G;
    difsums := [[k]];

    # perform search and remove equivalents for each stage
    for i in [2..Length(series)-1] do

        newN := series[i];
        difsums := SomeRefinedDifferenceSumsOptimized(G, oldN, newN, difsums);
        difsums := EquivalentFreeListOfDifferenceSums(G, newN, difsums);

        oldN := newN; # prepare for next stage

    od;

    # in final stage refine to difsets and remove equivalents
    # difsets := SomeRefinedDifferenceSets(G, oldN, difsums);
    difsets := SomeRefinedDifferenceSets(G, oldN, difsums);
    difsets := SmallestEquivalentFreeListOfDifferenceSets(G, difsets);

    return difsets;
end );


#############################################################################
##
#F  DifferenceSetsOfSizeKWithMultiplier( <G>, <k> )
##
InstallGlobalFunction( DifferenceSetsOfSizeKWithMultiplier, function (G, k)
    local series, oldN, difsums, i, newN, difsets;

    series := RefiningSeries(G);

    # at start the quotient group G/G is trivial and the only difsum is k
    oldN := G;
    difsums := [[k]];

    # perform search and remove equivalents for each stage
    for i in [2..Length(series)-1] do

        newN := series[i];
        difsums := SomeRefinedDifferenceSums(G, oldN, newN, difsums);
        difsums := EquivalentFreeListOfDifferenceSums(G, newN, difsums);

        oldN := newN; # prepare for next stage

    od;

    # in final stage refine to difsets and remove equivalents
    # difsets := SomeRefinedDifferenceSets(G, oldN, difsums);
    difsets := SomeRefinedDifferenceSetsWithMultiplier(G, oldN, difsums);
    difsets := SmallestEquivalentFreeListOfDifferenceSets(G, difsets);

    return difsets;
end );


#############################################################################
##
#F  DifferenceSetsOfSizeKUseMultiplier( <G>, <k> )
##
InstallGlobalFunction( DifferenceSetsOfSizeKUseMultiplier, function (G, k)
    local series, oldN, difsums, v, lambda, i, newN, difsets;

    series := RefiningSeries(G);

    # at start the quotient group G/G is trivial and the only difsum is k
    oldN := G;
    difsums := [[k]];
    v := Size(G);
    lambda := k*(k-1)/(v-1);

    # perform search and remove equivalents for each stage
    for i in [2..Length(series)-1] do

        newN := series[i];
        difsums := SomeRefinedDifferenceSumsUseMultiplier(G, k, lambda, oldN, newN, difsums);
        # difsums := EquivalentFreeListOfDifferenceSums(G, newN, difsums);

        oldN := newN; # prepare for next stage

    od;

    # in final stage refine to difsets and remove equivalents
    # difsets := SomeRefinedDifferenceSets(G, oldN, difsums);
    difsets := SomeRefinedDifferenceSetsUseMultiplier(G, oldN, difsums, k);
    difsets := SmallestEquivalentFreeListOfDifferenceSets(G, difsets);

    return difsets;
end );



#############################################################################
##
#F  DifferenceSetPairsOfSizeK( <G>, <k1>, <k2>, <e>, <lambda> )
# FEATURE DifferenceSetPairsOfSizeK(G, k1, k2, e, lambda)
InstallGlobalFunction( DifferenceSetPairsOfSizeK, function (G, k1, k2, e, lambda)
    local v, series, oldN, newN, difsums, i, difsetpairs;
    series := RefiningSeries(G);
    v := Size(G);
    oldN := G;
    difsums := [[[k1],[k2]]];
    # Print("\n\n\nDifferenceSetPairsOfSizeK>>>当前搜索大小为(", v, ", ", k1, ", ", k2, ", ", e, ", ", lambda, "):", "\n");
    # Print("DifferenceSetPairsOfSizeK>>>当前子群的阶为|N|=", v, ", 当前差集和为", difsums, "\n");
    for i in [2..Length(series)-1] do
        newN := series[i];
        difsums := SomeRefinedDifferenceSumsOfDifferenceSetPair(G, oldN, newN, difsums, k1, k2, e, lambda);
        # Print("DifferenceSetPairsOfSizeK>>>当前子群的阶为|N|=", Size(newN), ", 当前差集和为", difsums, "\n");
        difsums := EquivalentFreeListOfDifferenceSumsOfDifferenceSetPairs(G, newN, difsums);
        # Print("DifferenceSetPairsOfSizeK>>>去同构后所得差集和为", difsums, "\n");
        oldN := newN;  # prepare for next stage
    od;

    # in final stage refine to difsetpairs and remove equivalents
    difsetpairs := SomeRefinedDifferenceSetPairs(G, oldN, difsums, e, lambda);
    difsetpairs := EquivalentFreeListOfDifferenceSetPairs(G, difsetpairs);
    return difsetpairs;
end );


#############################################################################
##
#F  DifferenceSetPairsOfSizeKOptimized( <G>, <k1>, <k2>, <e>, <lambda> )
# FEATURE DifferenceSetPairsOfSizeKOptimized(G, k1, k2, e, lambda)
InstallGlobalFunction( DifferenceSetPairsOfSizeKOptimized, function (G, k1, k2, e, lambda)
    local v, series, oldN, newN, difsums, i, difsetpairs;
    series := RefiningSeries(G);
    v := Size(G);
    oldN := G;
    difsums := [[[k1],[k2]]];
    for i in [2..Length(series)-1] do
        newN := series[i];
        difsums := SomeRefinedDifferenceSumsOfDifferenceSetPairOptimized(G, oldN, newN, difsums, k1, k2, e, lambda);
        difsums := EquivalentFreeListOfDifferenceSumsOfDifferenceSetPairs(G, newN, difsums);
        oldN := newN;  # prepare for next stage
    od;

    # in final stage refine to difsetpairs and remove equivalents
    difsetpairs := SomeRefinedDifferenceSetPairsWithMultiplier(G, oldN, difsums, e, lambda);
    difsetpairs := EquivalentFreeListOfDifferenceSetPairs(G, difsetpairs);
    return difsetpairs;
end );

#############################################################################
##
#F  DifferenceSetPairsOfSizeKWithMultiplier( <G>, <k1>, <k2>, <e>, <lambda> )
# FEATURE DifferenceSetPairsOfSizeKWithMultiplier(G, k1, k2, e, lambda)
InstallGlobalFunction( DifferenceSetPairsOfSizeKWithMultiplier, function (G, k1, k2, e, lambda)
    local v, series, oldN, newN, difsums, i, difsetpairs;
    series := RefiningSeries(G);
    v := Size(G);
    oldN := G;
    difsums := [[[k1],[k2]]];
    # Print("\n\n\DifferenceSetPairsOfSizeKWithMultiplier>>>当前搜索大小为(", v, ", ", k1, ", ", k2, ", ", e, ", ", lambda, "):", "\n");
    # Print("DifferenceSetPairsOfSizeKWithMultiplier>>>当前子群的阶为|N|=", v, ", 当前差集和为", difsums, "\n");
    for i in [2..Length(series)-1] do
        newN := series[i];
        difsums := SomeRefinedDifferenceSumsOfDifferenceSetPair(G, oldN, newN, difsums, k1, k2, e, lambda);
        # Print("DifferenceSetPairsOfSizeKWithMultiplier>>>当前子群的阶为|N|=", Size(newN), ", 当前差集和为", difsums, "\n");
        difsums := EquivalentFreeListOfDifferenceSumsOfDifferenceSetPairs(G, newN, difsums);
        # Print("DifferenceSetPairsOfSizeKWithMultiplier>>>去同构后所得差集和为", difsums, "\n");
        oldN := newN;  # prepare for next stage
    od;

    # in final stage refine to difsetpairs and remove equivalents
    difsetpairs := SomeRefinedDifferenceSetPairsWithMultiplier(G, oldN, difsums, e, lambda);
    difsetpairs := EquivalentFreeListOfDifferenceSetPairs(G, difsetpairs);
    return difsetpairs;
end );

#############################################################################
##
#F  DivisibleDifferenceSetsOfSizeK( <G>, <w>, <v>, <l>, <lambda1>, <lambda2> )
# FEATURE DivisibleDifferenceSetsOfSizeK(G, w, v, l, lambda1, lambda2)
InstallGlobalFunction( DivisibleDifferenceSetsOfSizeK, function (G, w, v, l, lambda1, lambda2)
    local series, oldN, newN, difsums, i, diviDiffSets;
    series := RefiningSeries(G);
    oldN := G;
    difsums := [[l]];
    # Print("DivisibleDifferenceSetsOfSizeK>>>>当前搜索大小为(",w,", ",v,", ",l,", ",lambda1,", ",lambda2,")", "\n");
    # Print("DivisibleDifferenceSetsOfSizeK>>>>当前子群的阶为|N|=", w*v, ", 当前差集和为", difsums, "\n");
    for i in [2..Length(series)-1] do
        newN := series[i];
        difsums := SomeRefinedDifferenceSumsOfDivisibleDifferenceSet(G, oldN, newN, difsums, w, v, l, lambda1, lambda2);
        # Print("DivisibleDifferenceSetsOfSizeK>>>>当前子群的阶为|N|=", Size(newN), ", 当前差集和为", difsums, "\n");
        difsums := EquivalentFreeListOfDifferenceSumsOfDivisibleDifferenceSets(G, newN, difsums);
        # Print("DivisibleDifferenceSetsOfSizeK>>>>去同构后所得差集和为", difsums, "\n");
        oldN := newN;  # prepare for next stage
    od;

    # in final stage refine to diviDiffSets and remove equivalents
    diviDiffSets := SomeRefinedDivisibleDifferenceSets(G, oldN, difsums, w, v, l, lambda1, lambda2);
    diviDiffSets := EquivalentFreeListOfDivisibleDifferenceSets(G, diviDiffSets);
    return diviDiffSets;
end );

#############################################################################
##
#F  DivisibleDifferenceSetPairsOfSizeK( <G>, <w>, <v>, <l1>, <l2>, <le>, <lambda1>, <lambda2> )
# FEATURE DivisibleDifferenceSetPairsOfSizeK(G, w, v, l1, l2, le, lambda1, lambda2)
InstallGlobalFunction( DivisibleDifferenceSetPairsOfSizeK, function (G, w, v, l1, l2, le, lambda1, lambda2)
    local series, oldN, newN, difsums, i, diviDiffSetPairs;
    series := RefiningSeries(G);
    oldN := G;
    difsums := [[[l1],[l2]]];
    for i in [2..Length(series)-1] do
        newN := series[i];
        difsums := SomeRefinedDifferenceSumsOfDivisibleDifferenceSetPair(G, oldN, newN, difsums, w, v, l1, l2, le, lambda1, lambda2);
        difsums := EquivalentFreeListOfDifferenceSumsOfDivisibleDifferenceSetPairs(G, newN, difsums);
        oldN := newN;  # prepare for next stage
    od;

    # in final stage refine to diviDiffSetPairs and remove equivalents
    diviDiffSetPairs := SomeRefinedDivisibleDifferenceSetPairs(G, oldN, difsums, w, v, l1, l2, le, lambda1, lambda2);
    diviDiffSetPairs := EquivalentFreeListOfDivisibleDifferenceSetPairs(G, diviDiffSetPairs);
    return diviDiffSetPairs;
end );

#############################################################################
##
#F  DifferenceSets( <G> )
##
InstallGlobalFunction( DifferenceSets, function (G)
    local difsets, k; 

    difsets := []; # will store all inequivalent difsets

    for k in PossibleDifferenceSetSizes(G) do
        Append(difsets, DifferenceSetsOfSizeK(G, k));
    od;

    return difsets;
end );

#############################################################################
##
#F  DifferenceSetsOptimized( <G> )
##
InstallGlobalFunction( DifferenceSetsOptimized, function (G)
    local difsets, k; 

    difsets := []; # will store all inequivalent difsets

    for k in PossibleDifferenceSetSizes(G) do
        Append(difsets, DifferenceSetsOfSizeKOptimized(G, k));
    od;

    return difsets;
end );

#############################################################################
##
#F  DifferenceSetsWithMultiplier( <G> )
##
InstallGlobalFunction( DifferenceSetsWithMultiplier, function (G)
    local difsets, k; 

    difsets := []; # will store all inequivalent difsets

    for k in PossibleDifferenceSetSizes(G) do
        Append(difsets, DifferenceSetsOfSizeKWithMultiplier(G, k));
    od;

    return difsets;
end );

#############################################################################
##
#F  DifferenceSetsUseMultiplier( <G> )
##
InstallGlobalFunction( DifferenceSetsUseMultiplier, function (G)
    local difsets, k; 

    difsets := []; # will store all inequivalent difsets

    for k in PossibleDifferenceSetSizes(G) do
        Append(difsets, DifferenceSetsOfSizeKUseMultiplier(G, k));
    od;

    return difsets;
end );


#############################################################################
##
#F  DifferenceSetPairs( <G> )
##
# FEATURE DifferenceSetPairs(G)
InstallGlobalFunction( DifferenceSetPairs, function (G)
    local difsetpairs, k, dsp; 
    difsetpairs := []; # will store all inequivalent difsetpairs
    for k in PossibleDifferenceSetPairSizes(G) do
        dsp := DifferenceSetPairsOfSizeK(G, k[1], k[2], k[3], k[4]);
        Append(difsetpairs, dsp);
    od;
    return difsetpairs;
end );

#############################################################################
##
#F  DifferenceSetPairsOptimized( <G> )
##
# FEATURE DifferenceSetPairsOptimized(G)
InstallGlobalFunction( DifferenceSetPairsOptimized, function (G)
    local difsetpairs, k, dsp, order, totalCnt, f, map, i, size; 
    # difsetpairs := []; # will store all inequivalent difsetpairs
    totalCnt := 0;
    f := "";
    map := ["order","-","k1","-","k2","-","e","-","lambda",".txt"];
    order := String(Size(G));
    map[1] := order;
    for k in PossibleDifferenceSetPairSizes(G) do
        map[3] := String(k[1]);
        map[5] := String(k[2]);
        map[7] := String(k[3]);
        map[9] := String(k[4]);
        f := Concatenation(map);
        dsp := DifferenceSetPairsOfSizeKOptimized(G, k[1], k[2], k[3], k[4]);
        size := Size(dsp);
        # 搜完一个参数，就写入文件中
        for i in [1..size] do
            AppendTo(f, dsp[i]);
            if i< size then AppendTo(f, ","); fi;
            AppendTo(f, "\n");
        od;
        # Append(difsetpairs, dsp);
        totalCnt := totalCnt + size;
    od;
    return totalCnt;
end );

#############################################################################
##
#F  DifferenceSetPairsWithMultiplier( <G> )
##
# FEATURE DifferenceSetPairsWithMultiplier(G)
InstallGlobalFunction( DifferenceSetPairsWithMultiplier, function (G)
    local difsetpairs, k, dsp; 
    difsetpairs := []; # will store all inequivalent difsetpairs
    for k in PossibleDifferenceSetPairSizes(G) do
        dsp := DifferenceSetPairsOfSizeKWithMultiplier(G, k[1], k[2], k[3], k[4]);
        Append(difsetpairs, dsp);
    od;
    return difsetpairs;
end );

#############################################################################
##
#F  DifferenceSetPairsUseMultiplierOfAbelianGroupOfSizeK( G, k1, k2, e, lambda )
##
# FEATURE DifferenceSetPairsUseMultiplierOfAbelianGroupOfSizeK(G, k1, k2, e, lambda)
InstallGlobalFunction( DifferenceSetPairsUseMultiplierOfAbelianGroupOfSizeK, function (G, k1, k2, e, lambda)
    local difsetpairs, dif, v, multipliers,multiplier,tracksOfGroup, D1s, D1, D2s, D2; 
    difsetpairs := []; # will store all inequivalent difsetpairs
    dif := DifferenceTable(G);
    v := Size(G);
    multipliers := Multipliers(G, e, lambda);

    for multiplier in multipliers do
        tracksOfGroup := TracksOfGroup(G, multiplier);
        D1s:=SubsetSum(List([1..Size(tracksOfGroup)],x->Size(tracksOfGroup[x])), k1); # index
        D1s:=List(D1s, x->List(x, y->tracksOfGroup[y]));
        D1s:=List(D1s, x->Flat(x));

        D2s:=SubsetSum(List([1..Size(tracksOfGroup)],x->Size(tracksOfGroup[x])), k2); # index
        D2s:=List(D2s, x->List(x, y->tracksOfGroup[y]));
        D2s:=List(D2s, x->Flat(x));

        for D1 in D1s do
            for D2 in D2s do
                if Length(Intersection(D1,D2)) <> e then
                    continue;
                fi;

                if IsDifferenceSetPairByTable(dif, D1, D2, v, e, lambda) then
                    Add(difsetpairs, [D1,D2]);
                    break; # break when we found a dsp
                fi;
            od;
        od;

    od;
    difsetpairs := EquivalentFreeListOfDifferenceSetPairs(G, difsetpairs);
    return difsetpairs;
end );


#############################################################################
##
#F  DifferenceSetPairsUseMultiplierOfAbelianGroup( <G> )
##
# FEATURE DifferenceSetPairsUseMultiplierOfAbelianGroup(G)
InstallGlobalFunction( DifferenceSetPairsUseMultiplierOfAbelianGroup, function (G)
    local difsetpairs, dif, v, k, k1, k2, e, lambda, multipliers,multiplier,tracksOfGroup, D1s, D1, D2s, D2, dsp; 
    if not IsAbelian(G) then
        Print("This function is intended to search for Abelian Groups only!\n ");
        Print("Please provide an Abelian Group for searching!\n ");
        return [];
    fi;
    difsetpairs := []; # will store all inequivalent difsetpairs
    dif := DifferenceTable(G);
    v := Size(G);
    for k in PossibleDifferenceSetPairSizes(G) do
        k1 := k[1];
        k2 := k[2];
        e := k[3];
        lambda := k[4];
        multipliers := Multipliers(G, e, lambda);
        dsp := [];

        for multiplier in multipliers do
            tracksOfGroup := TracksOfGroup(G, multiplier);
            D1s:=SubsetSum(List([1..Size(tracksOfGroup)],x->Size(tracksOfGroup[x])), k1); # index
            D1s:=List(D1s, x->List(x, y->tracksOfGroup[y]));
            D1s:=List(D1s, x->Flat(x));

            D2s:=SubsetSum(List([1..Size(tracksOfGroup)],x->Size(tracksOfGroup[x])), k2); # index
            D2s:=List(D2s, x->List(x, y->tracksOfGroup[y]));
            D2s:=List(D2s, x->Flat(x));

            for D1 in D1s do
                for D2 in D2s do
                    if Length(Intersection(D1,D2)) <> e then
                        continue;
                    fi;

                    if IsDifferenceSetPairByTable(dif, D1, D2, v, e, lambda) then
                        Add(dsp, [D1,D2]);
                        break; # break when we found a dsp
                    fi;
                od;
            od;

        od;
        dsp := EquivalentFreeListOfDifferenceSetPairs(G, dsp);
        Append(difsetpairs, dsp);
    od;
    return difsetpairs;
end );

#############################################################################
##
#F  DivisibleDifferenceSets( <G> )
##
# FEATURE DivisibleDifferenceSets(G)
InstallGlobalFunction( DivisibleDifferenceSets, function (G)
    local dividifsets, k, dds; 
    dividifsets := [];
    for k in PossibleDivisibleDifferenceSetSizes(G) do
        dds := DivisibleDifferenceSetsOfSizeK(G, k[1], k[2], k[3], k[4], k[5]);
        Append(dividifsets, dds);
    od;
    return dividifsets;
end );


#############################################################################
##
#F  DivisibleDifferenceSetPairs( <G> )
##
# FEATURE DivisibleDifferenceSetPairs(G)
InstallGlobalFunction( DivisibleDifferenceSetPairs, function (G)
    local dividifsetpairs, k, ddsp; 
    dividifsetpairs := [];
    for k in PossibleDivisibleDifferenceSetPairSizes(G) do
        ddsp := DivisibleDifferenceSetPairsOfSizeK(G, k[1], k[2], k[3], k[4], k[5], k[6], k[7]);
        Append(dividifsetpairs, ddsp);
    od;
    return dividifsetpairs;
end );

#############################################################################
##
#F  DifferenceSumsOfSizeK( <G>, <N>, <k> )
##
InstallGlobalFunction( DifferenceSumsOfSizeK, function (G, N, k)
    local phi, series, oldN, difsums, i, newN;

    # produce a refining series for G/N and find its preimage in G
    phi := NaturalHomomorphismByNormalSubgroup(G, N);
    series := RefiningSeries(Image(phi));
    Apply(series, x -> PreImagesSet(phi, x));

    # at start the quotient group G/G is trivial and the only difsum is k
    oldN := G;
    difsums := [[k]];

    # perform search and remove equivalents for each stage
    for i in [2..Length(series)] do

        newN := series[i];
        difsums := SomeRefinedDifferenceSums(G, oldN, newN, difsums);
        difsums := EquivalentFreeListOfDifferenceSums(G, newN, difsums);

        oldN := newN; # prepare for next stage

    od;

    return difsums;
end );

#############################################################################
##
#F  DifferenceSums( <G>, <N> )
##
InstallGlobalFunction( DifferenceSums, function (G, N)
    local difsums, k;

    difsums := []; # will store all inequivalent difsums

    for k in PossibleDifferenceSetSizes(G) do
        Append(difsums, DifferenceSumsOfSizeK(G, N, k));
    od;

    return difsums;
end );

#############################################################################
##
#E

